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Charge transport through a nanoscale junction coupled to two macroscopic electrodes is investi- 
gated for the situation when bound states are present. We provide numerical evidence that bound 
states give rise to persistent, non-decaying current oscillations in the junction. We also show that 
the amplitude of these oscillations can exhibit a strong dependence on the history of the applied 
potential as well as on the initial equilibrium configuration. Our simulations allow for a quantitative 
investigation of several transient features. We also discuss the existence of different time-scales and 
address their microscopic origin. 
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I. INTRODUCTION 

In order to describe electronic transport through meso- 
scopic or nanoscopic devices, a quantum description of 
transport is essential. A seminal quantum theory of 
transport is the Landauer-Biittiker formalism,— » which 
expresses the conductance of a device in terms of the 
quantum-mechanical transmittance of (non-interacting) 
electrons at the Fermi energy. 

In recent years and spurred by experimental 
progress in transport measurements through single 
molecules, 3 the Landauer-Biittiker formalism has been 
combined^'^'^^ ' 10 ! 11 ] 12 with (static) density functional 
theory which allows to take the atomistic structure of 
both the molecule and the contacts into account. For a 
recent critical review of this methodology, the reader is 
referred to Ref. 

The Landauer-Biittiker formalism focusses on the de- 
scription of steady-state transport and assumes that for a 
system which is driven out of equilibrium by a dc bias, a 
dc current will eventually develop, which means that the 
dynamical formation of the steady state is not proved but 
rather taken for granted. The question how the system 
dynamically reaches a steady state has been investigated 
both numericall y 14 ' 15 ! 16 ! 17 and theoretically 18 ' 19 Using 
non-equilibrium Green functions (NEGF) techniques it 
has been shown-?, that the total current (and density) 
approaches a steady value provided the local density of 
states is smooth in the device region. Such value is 1) 
in agreement with the Landauer formula and 2) indepen- 
dent of the initial equilibrium configuration and the his- 
tory of the applied bias. For a steady state to develop the 
condition on the local density of states excludes the pres- 
ence of bound states. Recently, the inclusion of bound 
states in time-dependent quantum transport has been 
studied in Ref. [2fJ and further been addressed in sub- 
sequent workiSL There it is demonstrated that if the dc 
biased Hamiltonian supports two or more bound states, 
the long-time limit of the current consists of two terms: 
a steady-state contribution given by the Landauer for- 



mula and an additional, dynamical contribution respon- 
sible for undamped current oscillations. The frequencies 
of these oscillations are given by the differences between 
two bound-state energies and, interestingly, the ampli- 
tudes depend on both the initial state and history of the 
time-dependent perturbation. 

In the present work, the history as well as the initial- 
state dependence of the dynamical part of the current 
is investigated numerically in detail. As a tool for 
our numerical calculations we use a recently developed 
algorithm 1 ^ which allows for the time propagation of 
quantum transport systems according to the Schrodinger 
equation. 

The paper is organized as follows. In Section [TTJ we 
summarize the results of Ref. [2l| which are relevant for 
the discussion of our findings and we briefly describe the 
central ideas of the time-propagation algorithm. In Sec- 
tion IIIII we present our numerical results which not only 
confirm the existence of the undamped current oscilla- 
tions but also allow to identify additional internal transi- 
tions contributing to the transient behavior of the driven 
system. We investigate the dependence of the current 
oscillations on various parameters and initial conditions 
and provide theoretical explanations of the observed be- 
havior. Finally, we recapitulate our main results in Sec- 
tion us 



II. TWO APPROACHES TO 
TIME-DEPENDENT TRANSPORT 

In this Section we briefly describe two alternative 
approaches to time-dependent transport in a typi- 
cal electrode-device-electrode geometry: non-equilibrium 
Green functions (NEGF) and direct solution of the time- 
dependent Schrodinger equation. As was already pointed 
out within the former approach ; 20 ' 21 quantum transport 
in systems of non-interacting electrons exhibits persis- 
tent current- (and density-) oscillations if two or more 
bound states are present in the biased system. Here, we 
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use the latter approach to address several issues about 
such bound-state oscillations. A particularly interesting 
feature of them is the fact that their amplitude depends 
on the entire time evolution as the system is driven out 
of equilibrium (memory effects). 



A. Non-equilibrium Green functions 

We consider a quantum system of non-interacting elec- 
trons which consists of a central device (e.g., a quantum 
point contact or a single molecule plus a few atomic lay- 
ers of the left and right electrodes) and two semi-infinite 
reservoirs (left and right electrodes). As initial state we 
use the one proposed by Cinij^ all parts of the system, 
i.e., left lead (region L), central device (region C) and 
right lead (region R), are initially (at t < 0) connected 
and in a well defined equilibrium configuration with a 
unique temperature and chemical potential (thermody- 
namic consistency). In this initial state, the charge den- 
sity of the electrodes is perfectly balanced and no current 
flows through the junction. 

For non-interacting electrons at zero temperature, the 
initial state is a Slater determinant of eigenstates of the 
entire contacted system with eigenenergies smaller than 
the Fermi energy. At time t > the system is driven 
out of equilibrium by exposing it to an external time- 
dependent potential which is local in time and space. 
For example, we may switch on an electromotive force in 
such a way that the potential drop is entirely limited to 
the central region. The boundaries of the open quantum 
system are chosen in a way that the density outside the 
region C is accurately described by an equilibrium bulk 
density. The time-dependent perturbation may cause a 
current flow through the device. The total current from 
region a = L, R can be calculated from time derivative 
of the total number of particles in a: 



I a (t) 



L, R, 



(1) 



where n(r,t) is the time-dependent electron density and 
the space integral extends over region a (e is the electron 
charge). Assuming no direct coupling between the left 
and right electrodes, the single-particle Hamiltonian of 
the entire, contacted system can be written as: 



H(t) 



H LL (t) H LC 
Hcl H C c(t) H CR 
H RC H RR (t) 



(2) 



The diagonal blocks of the above matrix are obtained by 
projecting the full Hamiltonian H onto the correspond- 
ing region. The off-diagonal blocks in Eq. @ account 
for the coupling between the device region C and the 
leads and, for simplicity, we assume them to be time- 
independent. For instance, in a real-space representation 
using a finite-difference discretization of the kinetic en- 
ergy, the off-diagonal elements of H are simply given 



by the off-diagonal elements of the kinetic energy opera- 
tor. (Model systems with time-dependent couplings were 
studied, e.g., in Ref. I22I) 

One way to deal with non-equilibrium problems is pro- 
vided by the NEGF theory. From the equation of motion 
of the Keldysh-Green function one can rewrite the cur- 
rent I a (t) of Eq. |T]) in terms of the lesser Green function 
projected onto different subregions as: 



I a (t) = 2eReTr[G% a (t,t)H aC ], 



(3) 



where Tr denotes the trace over a complete set of states 
in the central region. The lesser Green function can be 
expresse d 18 i 19 ' 23 i 24 ' 25 in terms of retarded and advanced 
Green functions as 



G < (t;t') = G R (t;0)G < (0;0)G A (0;t'). 



(4) 



The initial condition is G<(0; 0) = if(H°) where f(u) = 
(e^^^ + l) -1 is the Fermi distribution function and H° 
is the (time-independent) Hamiltonian for t < 0. 

It can be shown 2 ^ that in a dc biased system the total 
time-dependent current approaches a steady value pro- 
vided the local density of states in region C is smooth. In 
this case, the steady current is given by: 

l[ S) = l}ml a (t) = e J ^\f( u -U?)-f(w-U!t)]T(u). 

(5) 

In the above equation is the value ap- 

proached by the bias in lead a when t — > 00 and 
T(u) = Tt[G^ c (u)T L (u)G^ c (u)T R (u)], where 
T a (uj) — ~ 2Im[Sf (uj)] and G R ^ are the retarded 
and advanced Green functions projected in region 
C. Sf (uj) = Hc a g R a {u)H a c is the embedding self 
energy with the retarded Green function of lead a, 



Hi 



U c 



i0+) 



The steady 



current does not depend on the initial Hamiltonian 
(the memory of different initial conditions is completely 
washed out) and is also independent of the history of 
the applied bias (memory-loss theorem)?^ 

The above scenario changes drastically if the Hamilto- 
nian H°° := lim^oo H{t) has two or more bound eigen- 
states. In this case the long-time limit of the current has 
two contributions: 21 



lim I a (t) = 1^ + li D \t) 



(6) 



(S) 

In addition to the steady-state contribution la given by 
Eq. ([5]) one finds a dynamical, explicitly time-dependent 
contribution I a D ^ which can be written as 



')*]■ (7) 



6,6' 



In Eq. ([7]) the summation is over all bound states of the fi- 
nal Hamiltonian H 00 and I a D ' oscillates with frequencies 
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given by the differences of the bound-state eigenenergies. 
The quantities Abb' an d fb,b' are defined according to 



A 



b,b> 



Tr r 



m)Wc\vi{t' b ) 



and 



h,v = mf{H°)w b ,) 



(8) 



(9) 



The state iV'fcc) is the projection of the bound eigenstate 
iV'ft ) of the biased Hamiltonian H°° onto the central 
region. The state \ip' b ) is related to IV'? ) by a unitary 
transformation: 



W bL ) 




~e lA ?l L 






W bC ) 




M c 






Wm) 




e iA *l R 







with 



A~ = lim 

t— »oo 



dt'(U a (f ) - t/~), 



(10) 



(11) 



Mc a unitary "memory matrix" with the same dimen- 
sion as the number of degrees of freedom employed to 
describe region C and l a the identity matrix projected 
onto region a — L, R . The memory matrix depends 
on the history of the time-dependent perturbation and is 
defined through the equation below 



lim G£ c (0;t) = M c lim G cc (0;t), 

t — >oo t — >oo 



(12) 



where Gp C (0;i) is the projection onto region C of the 
advanced Green function G A (0;i) = iexp(iH°°t). 

Few remarks about the central result in Eq. are in 
order. First, we wish to emphasize again that no steady- 
state current develops if the biased Hamiltonian H°° 
has bound eigenstates. The current oscillations given by 
Eq. ([7]) are persistent, i.e., they do not decay in time. 
Second, in contrast to the case without bound states, the 
asymptotic current I a (i) depends both on the initial equi- 
librium configuration and history of the applied bias and 
gate voltage through the coefficients fb,v of Eq. ©. For 
sudden switching of the bias and gate voltage A^° = 
and Mc = lc (lc being the identity matrix projected 
onto region C) and the matrix in Eq. (fit)]) reduces to 
the identity matrix. On the contrary, different switch- 
ing processes yield different memory matrices and hence 
different amplitudes of the current oscillations, see Sec- 
tion IIIII for a detailed study of the history dependence. 
Third, the NEGF formalism described in this Section can 
be combined with Time-Dependent Density Functional 
Theory2&2I (TDDFT) to include exchange and correla- 
tion effects in the calculated density and current. In this 
theory the steady-state assumption is consistent with the 
TDDFT equation for the total current provided the den- 
sity of states in region C is a smooth function^ On 
the contrary, the presence of bound-states in the biased 
Hamiltonian is not compatible with a steady current 3i 



This result opens up the possibility of having oscillatory 
solutions even for constant biases and may change sub- 
stantially the standard steady-state picture already at 
the level of exchange-correlation functionals which are 
local or semi-local in time. On one hand, the oscilla- 
tions of the effective potential in region C give rise to 
new conductive channels, an effect that cannot be cap- 
tured in any static approach. On the other hand, the 
asymptotic (t — > oo) density depends on the occupation 
coefficients fb y b which in turn depend on the history of 
the TDDFT potential. Thus, history-dependent effects 
might be observed even at the level of the adiabatic lo- 
cal density approximation. Finally we emphasize that 
the above conclusions are not limited to TDDFT but 
also apply to any other single-particle theory of electrons 
such as, e.g., Hartree-Fock theory. Similarly, they also 
apply to a single-electron theory of coupled electronic 
and nuclear motion where the time evolution of the nu- 
clei is treated in the Ehrenfest approximation and thus 
the potential acting on the electrons depends paramet- 
rically on the (time-dependent) nuclear coordinates. In 
this latter case the presence of a self-consistent oscillatory 
solution in a Holstein wire connected to one-dimensional 
non-interacting leads was observed in Ref. |28| . 



B. Direct propagation of the time-dependent 
Schrodinger equation 

Calculating the time-dependent current in terms of 
the Green function projected onto the central region 
amounts to solving either the Keldysh-Dyson integral 



equation, 
equations 



,29,30 
31,32 



or the integro-differential Kadanoff-Baym 
In this work we use an alternative ap- 
proach which is based on solving the time-dependent 
Schrodinger equation for the initially occupied one- 
particle statesJ^ An advantage of the latter approach 
over the former ones is that the wave-functions depend 
only on one time argument as opposed to the double time 
dependence of the Green function. This algorithm has re- 
cently also been used to study electron pumping by direct 
time propagation^. 

For non-interacting electrons at zero temperature the 
total current from region a of Eq. (JT|) can alternatively 
be expressed as a surface integral 



Ia(t) 



/ dan-lm[>P*(r,t)VTp n (r,t)}, (13) 

J Sr. 



where n is the unit vector perpendicular to the surface el- 
ement da, the surface S a is perpendicular to the longitu- 
dinal geometry of our system and ip n (r, 0) are the eigen- 
states of H(t < 0). The electrode-junction-electrode sys- 
tem is infinitely extended and non-periodic. In practice, 
of course, we can only deal with finite systems and there- 
fore we only propagate the initial wavefuction projected 
onto the central region C. The presence of the leads 
is taken into account by applying the correct boundary 
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conditions. It is worth to note that even for interact- 
ing electrons one can use Eq. (fT5|) to compute the cur- 
rent through the junction if the single-particle orbitals 
ip n (r,t) are the Kohn-Sham orbitals of time-dependent 
density functional theory. 

For a description of the algorithm proposed in Ref. [Til , 
it is convenient to write H aa (t), with a = L,R , as 
the sum of a term H a aa — H aa (0) which is constant in 
time and another term U a (t) which may be explicitly 
time-dependent, H aa (t) = H aa + U a (t). In configu- 
ration space U a (t) is diagonal at any time t since the 
potential is local in space. Furthermore, the diagonal 
elements U a (r, t) are spatially constant for metallic elec- 
trodes. Thus, U a {t) = U a (t)l a and U L (t) - U R (t) is 
the total potential drop across the junction. The total 
Hamiltonian is H(t) = H{t) + U(t) with 



H(t) 



and 



U(t) = 



Hcl Hcc(t) Hen. 
H RC H rr 



U L {t)l L 

U R (t)l R 



(14) 



| ( S , v m )) depends on the initial wavefunction in region 
a = L, R and reads 



2iS 



iSH 



E 



(to) 

eff a=L, R 



»(m,0) 



H 



Ca 



(i-iSH aa y 

(l + iSH aa ) m 



with 



,(«>) - 



A S Jj( m ) 

- " and kt> k) 



i+4ui m) 



j=k 



(0)\ 



,0-)l 2 



(17) 



(18) 



The memory term \M^ m ^} is responsible for the hopping 
in and out of region C . It depends on the wavefunction 
in the device region at previous time steps and reads 



TO— 1 i (m,k) 



E 



1 + wil cS a=L, R k=0 Ua Ua 



(m — k) 



Q 



(m—k 



(\4 k+1) ) + l^>) . (19) 



withQi m) = H Ca [(l-i8H aa ) m /{l + i8H aa ) m + 1 ]H aC . 
For more details on the implementation of the algorithm 
the reader is referred to Ref. uA, 



In this way, the only term in H(t) that depends on 
t is Hcc(t)- For any given initial one-particle state 
|,/;(0)) = |^(°)> we calculate \i>(t m = mAt)) = \^) 
by employing a generalized form of the Cayley method 
(atomic units are used throughout) 



1 + iSH 



("i 



i + 4c/ (m) 



IV' 



(m+l)\ _ 



(m) 



(to)\ 



(15) 



with H ■* 



i[H(t m+1 )+ff(t m )], = \[U(t m+1 ) + 
U(t m )] and S = At/2. The above propagation scheme is 
unitary (norm conserving) and accurate to second-order 
in 6. From Eq. (|15jl we can extract an equation for the 
time-evolved state in region C . After some algebra, one 
ends up with an equation which gives the wave function 
in region C at time step m + 1 in terms of the wave 
function in region C at the previous time step and two 
additional terms (source and memory term): 



I , (m+lV 



la - iSH 



(to) 



lc + iSH 



^|4 m )) + |<?M)-|MM). (16) 



The effective Hamiltonian of region C is de- 

fined according to = Hqq — i5Hcl{^ + 

i&H\ L y l H hC - iSHcnO- + iSH^^Hac, where 
H CC = h[ H cc{Un+i) + H C c(t m )]. The source term 



III. NUMERICAL RESULTS 

In this Section we present the results of our numer- 
ical simulations for simple one-dimensional model sys- 
tems which support two bound states in the long-time 
limit. Of particular interest will be the dynamical part 
of the current and the dependence of the amplitude of 
the bound-state oscillations on the history of the time- 
dependent potential and on the initial state. We also 
identify single-particle transitions other than between the 
bound states which are relevant to understand the shape 
of the transient current. 

The time-dependent, one-dimensional Hamiltonian is 
given by 

1 d 2 

H(x,t) = --—s+Uo(x) + U(x,t) =: H°(x) + U(x,t) . 

2 dar 

(20) 

For times t < the Hamiltonian is H° (x) and the system 
is in its ground state. At t = the system is driven out of 
equilibrium by the time-dependent potential U(x, t). We 
choose the time-dependent perturbation in such a way 
that for t — > oo the Hamiltonian globally converges to an 
asymptotic Hamiltonian, which we denote with H°°(x). 

The time-dependent perturbation U(x,t) can be writ- 
ten as a piece- wise function of the space variable x. Let 
U a (t) be the applied bias in region a = L,R and V g (x,t) 
the gate voltage applied to region C. The latter may 
depend on both position x and time t. Then 



UL(t) — oo < x < xl 
U(x,t)=< V g (x,t) X L < x < X R 
U R {t) x R < x < oo 



(21) 
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with xl and xr the positions of the left and right inter- 
faces respectively. In our numerical implementation we 
discretize H on a equidistant grid and use a simple three- 
point discretization for the kinetic energy. In all systems 
studied below the simulations have been performed by 
considering a propagation window which extends from 
Xl = —1.2 a.u. to xr — 1.2 a.u. and a lattice spacing 
Ax — 0.012 a.u.. The occupied part of the continuous 
spectrum ranges from k — to kp — \/2ep and it is 
discretized with 200 /c-points. All occupied states are 
propagated from t = to t = 1400 a.u. using a time step 
25 = 0.05 a.u.. In all the numerical examples studied 
below the final Hamiltonian supports two bound states 
and the resulting current in the long-time limit then is 

I(t) +I osc (x)sm(uj t) (22) 

and, on top of the steady current I^ s \ has an oscillating 
part with only one frequency ojq given by the eigenen- 
ergy difference of the two bound states. It is also worth 
mentioning that the amplitude 7 sc of this current oscil- 
lation depends on the position (see Eq. ([5])) while the 
steady-state current is position-independent. 

A. Bound state oscillations and transients 

As a first example, we study a system with an ini- 
tial potential Uq(x) — 0. Initially, the system is in the 
ground state with Fermi energy ep =0.1 a.u.. All wave- 
functions of the ground-state Slater determinant are ex- 
tended one-particle states with energy between and Ep. 
At t = 0, the system is suddenly driven out of equilib- 
rium by switching on a potential U(x,t) which consists 
of a constant bias in the left lead, Ul = 0.1 a.u., and a 
constant gate voltage in the central region, V g = —1.4 
a.u.. The biased Hamiltonian has two bound eigenstates 
with energies = —1.032 a.u. and eg° 2 = —0.133 a.u.. 
From the discussion of the previous Section we expect 
that a steady state cannot develop and that the time- 
dependent current exhibits an oscillatory behavior with 
frequency ujq = e£° 2 — ej*^. This is indeed confirmed by 
our numerical simulations, as one can see in Fig. [TJ where 
we plot the modulus of the discrete Fourier transform of 
the time-dependent current. The latter quantity is de- 
fined according to 

n p +N 

I(ui k ) = £ I(2n6)e-^, ui k = f* 

7TV2-/Vo ^ N 
n—n p 

We have computed I(uJk) for different values of n p = 
(4 + 2p) ■ 10 3 , p = 0, 1, 2, 3, 4, and N = 16 • 10 3 . Dif- 
ferent values of p correspond to different time intervals 
te (tp,t p + T ) with t p = (2 +p) x 100 a.u. but with the 
same duration To = 800 a.u.. The coefficient in Eq. ([23]) 
is defined such that the height of the peak I(ui) at ui is 
equal to the amplitude of the oscillations with frequency 
ui. Besides the zero- frequency peak (not shown) due to 



0.04 




<o /a.u. 

FIG. 1: Modulus of the discrete Fourier transform of the cur- 
rent for V g = —1.4 a.u. and a constant bias in the left lead 
Ul = 0.1 a.u.. The inset shows a magnification of the region 
with bound-continuum transitions from the bound state with 
higher energy to the Fermi energy. Different curves corre- 
spond to different time intervals. 



the non vanishing dc current, I(ui) shows a dominant 
peak at the frequency uiq — e£° 2 — of the transition 
between the two bound states. As expected, the height 
of this peak remains unchanged as p varies from to 4, 
i.e., the current oscillation associated with this transi- 
tion remains undamped. We emphasize that they are an 
intrinsic property of the biased system. 

Closer examination of Fig. [T] reveals four extra peaks 
which are related to different internal transitions. The 
first and the last pairs of peaks occur at frequencies which 
correspond to transitions between the bound states and 
the lower edge of the unoccupied part of the continuous 
spectrum in the left and right lead of the biased system, 
e 6°i — * £ F' an d e b°i "~ * £ f + Ul, with i = 1,2. These sharp 
structures (mathematically stemming from the disconti- 
nuity of the zero-temperature Fermi distribution func- 
tion) give rise to long-lived oscillations of the total cur- 
rent and density. These oscillatory transients die off very 
slowly, the height of the peaks decreases with increasing 
t p empirically as l/t p (power-law behavior). In Fig. [TJ 
as well as in all following examples, we report results for 
the current calculated in the center of the device region. 
However it is worth to mention that the amplitude of 
the current oscillations decays exponentially in the leads 
as e -(*r.i+*?, a )k-*«l where k^ = ^2(\e^\ + U a ) with 

i = 1, 2, a = L, R and a; is a point in lead a. Conse- 
quently, the dynamical part of the current vanishes deep 
inside the leads (away from where the bound states are 
localized). 

In the second example, we consider a system described 
by the translationally invariant Hamiltonian H(x,t < 
0) = — g jra. At t — we suddenly switch on a con- 
stant bias in the left lead Ul = 0.15 a.u. and propagate 
until T = 150 a.u. when a steady state is reached. At 
t = T a gate voltage V g (x) = —v g = —1.02 a.u. is 
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FIG. 2: Modulus of the discrete Fourier transform of the cur- 
rent of a translationally invariant initial Hamiltonian which is 
perturbed at t = by a sudden bias in the left lead Ul =0.15 
a.u. and the system evolves toward a steady state. Then, at 
T — 150 a.u. a gate voltage V g (x) = —v g — —1.02 a.u. is sud- 
denly turned on. The first peak appears at the u — 0.686 a.u. 
which is the modulus of the energy of the bound eigenstate of 
the final Hamiltonian (H(x, t > T) has one bound eigenstate). 
Different curves correspond to different Fermi energies. 



suddenly turned on and the Hamiltonian H(x,t > T) 
has one bound eigenstate at energy = —0.686 a.u. . 
The depth v g is chosen in such a way that if one slightly 
increases v g a second bound eigenstate appears. Since 
the system has only one bound state, the oscillations die 
out slowly as 1/ (t — T) and eventually another steady 
state develops. In order to understand the transient os- 
cillations we have studied the Fourier transform of the 
current as shown in Fig. [2l There the first peak appears 
at the frequency of u> = (e^l which is a transition be- 
tween the bound level and the bottom of the continuum. 
As such, the position of this peak remains unchanged 
for different Fermi energies. Besides this transition one 
observes other peaks whose positions shift as the Fermi 
energy is changed. They correspond to transitions from 
the bound level to the top of the left and right continua 
and, as for the first transition, they decay as l/(t — T). 



B. Dependence of the current oscillations on the 
initial conditions 

The dynamical part of the current depends on the ini- 
tial Hamiltonian H°(x) through the amplitudes fb.v of 
Eq. In the first example of the previous Section 

the Hamiltonian at negative times, H (x), had no bound 
eigenstates. At positive times a gate voltage and a bias 
in the left lead were suddenly switched on and the Hamil- 
tonian at positive times is equal to H°°(x) and has two 
bound eigenstates. We now consider a system with two 
bound eigenstates for t < and exposed to a dc bias 
for t > 0. Specifically, we start with a static poten- 
tial describing a quantum well of depth Uq(x) = —1.4 



a.u. for |a;| < 1.2 a.u.. The ground state of the sys- 
tem is the Slater determinant of all the extended eigen- 
states with energy up to £f = 0.1 a.u. and of the two 
bound eigenstates at energies 1 = —1.035 a.u. and 
e|j 2 = —0.156 a.u.. At t = a dc bias Ur = 0.1 a.u. is 
suddenly switched on in the left lead and the Hamilto- 
nian H(x,t > 0) = H°°(x) is equal to the final Hamil- 
tonian studied in the previous Section. The resulting 
time-dependent current for these two systems are shown 
in Fig. H 




100 200 

t/a.u. 



FIG. 3: Comparison of the time-dependent current for sys- 
tems with and without bound states at negative times. The 
inset shows a magnification of the time-dependent current of 
the system with two initial bound states. Since both sys- 
tems have the same final Hamiltonian, the frequencies of the 
current oscillations are the same while the amplitude of the 
oscillations for the quantum well (with two bound state ini- 
tially) is smaller by almost two orders of magnitude than for 
the system without initial bound states. 

As a consequence of the fact that H°° (x) is the same in 
both systems the time-dependent currents should oscil- 
late with the same frequency, a result which is confirmed 
by our numerical calculation. The amplitude of this oscil- 
lation, however, depends on the initial equilibrium config- 
uration as well as on how H(x, t) approaches the asymp- 
totic Hamiltonian H ca {x). As one can see from Fig. |31 
the amplitude is much larger in the system with no initial 
bound states. This difference can be explained qualita- 
tively by looking at Eq. In both systems the time- 
dependent perturbation is switched on suddenly. There- 
fore, the transformation matrix of Eq. (flU|) becomes the 
unit matrix and Eq. ([9]) reduces to 

f h # = mf(H°)\^). (24) 

When the perturbation is small like in the case of the 
system with two initial bound states (H° H°°), the 
eigenf unctions l^? ) of H°° are approximate cigenfunc- 
tions of H° as well. Therefore /(IT )|VO « /(£b)|i/> 6 °°) 
and fb.b' ~ f(£b)$b,b> which leads to a vanishing dynam- 
ical part of the current since there only the off-diagonal 
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FIG. 4: The amplitude of the current oscillation as function 
of the switching time of the bias. The bias in the left lead is 
switched according to Ul (t) = Ul sin 2 (w&i) for t < t b = ^j- 
and ?7l(4) = Ul =0.1 a.u. for later times. The frequency of 
the current oscillation wo = e^2 — e 6°i is given by the difference 
of bound state energies in the final system which have the 
values e^i = —0.933 a.u. and e^ 2 = —0.063 a.u., respectively. 
The Fermi energy is ep = 0.2 a.u. and the gate potential is 



Va 



-1.3 a.u. 



elements contribute. By contrast, if the applied poten- 
tial U(x,t) is large, the overlap (ip^\f(H°)\ip™) can be 
quite substantial and the resulting amplitude of the cur- 
rent oscillation is large. 



creases monotonically as a function of if,, a trend which 
is expected in the region of long switching times (adi- 
abatic switching). Such behavior, however, does not 
contradict the analytic results of Section III Al In fact, 
the memory matrix in the central region Mq also de- 
pends on the way the bias is switched on through the 
time-dependent embedding self-energy needed to calcu- 
T2"|) . and, in general, Mp ^ lc 

2ir, 47T, . . . 



late Gg c (0;i), see Eq 



when A|° 



D. Dependence of the current oscillations on the 
history of the gate voltage 

Finally we present some results to illustrate the depen- 
dence of the current oscillations on the switching process 
of the gate voltage. Again we start with the constant 



t<0 



t>T + t„ 



U, 





i 












/ 


e b,2 





AX =2.4 



1 



C. Dependence of the current oscillations on the 
history of the bias 

The amplitude of the bound state oscillations depends, 
through the transformation matrix in Eq. (|10p . on the 
history of the time-dependent potential which perturbs 
the initial state. In this Section we investigate for the 
first time how such amplitudes depend on the switching 
process (history-dependence effects). 

We take the flat potential Uq(x) — as initial potential 
and the Fermi energy cf = 0.2 a.u.. At t — a gate volt- 
age Vg(x) = — 1.3 a.u. abruptly lowers the potential in 
the center. In addition, a time-dependent bias is applied 
to the left lead as f/i,(i) = Ul sin 2 (w;,i) for i < if, = 
and U L (t) — U L for t > where U L = 0.1 a.u.. 

The final biased Hamiltonian has two bound states 
with energies e?j = —0.933 a.u. and eg° 2 = —0.063 a.u. 
which again leads to undamped oscillations in the cur- 
rent. 

Choosing if, in such a way that A|° equals 2ir, 
Ait,. . . the upper block of the unitary matrix in Eq. (|10p 
become the identity matrix in region L. This suggests 
that the amplitude of the current oscillations might ex- 
hibit a non-monotonic behavior as a function of the 
switching time. Our numerical results demonstrate that 
this is not the case. Fig. [J] shows that the amplitude de- 



FIG. 5: Schematic sketch of the time evolution of the Hamil- 
tonian. Starting from an initially constant potential (left), at 
t = a bias is suddenly applied to the left lead and the 
system evolves toward a steady state (center). Then, be- 
tween times T and T + t g , a time-dependent gate voltage 
V g {x, t) — — -f-(t — T) is switched on in region C. For times 
t > T + t g (right) the Hamiltonian remains constant in time. 

potential Uq(x) = at equilibrium. At t = a bias 
is ramped up abruptly in the left lead and the time- 
dependent current goes through some transient which 
lasts for a few tens of atomic units. We wait long enough, 
a time T = 150 a.u., for a steady-state to develop. After 
this time all dependence on the history of the applied 
bias is washed out. 

At i = T a time dependent gate voltage V g (x,t) = 
— T&Ct — T) is applied to region C. The gate voltage de- 
creases linearly until i = T + t g and remains constant 
and equal to — v g for all later times. In Fig. [5]we provide 
a schematic sketch of the overall time-dependent pertur- 
bation. 

The time t g is the switching time. The final Hamilto- 
nian H co (x) = H(x,t > T + t g ) has two bound eigen- 
functions and the steady-state cannot develop. 

In Fig. [5] the amplitude of the oscillation versus the 
switching time t g is shown for a final depth of the gate 
v g = 1.3 a.u.. In the upper panel, the bias in the left lead 
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t /a.u. 




FIG. 6: The amplitude of the current oscillations as function 
of the switching time t g for v g = 1.3 a.u.. Upper panel: for 
fixed bias Ul =0.15 a.u. and different Fermi energies. Lower 
panel: for fixed Fermi energy ef = 0.2 a.u. and different 
values of the bias. All curves reach a maximum whose position 
remains almost unchanged. 



is fixed to Ul = 0.15 a.u. and the Fermi energy is varied 
from £p = 0.1 a.u. to 0.3 a.u.. We see that the amplitude 
reaches a maximum value for a certain switching time. 
It is also worth noting that the amplitudes are generally 
smaller for larger Fermi energies, a behavior which can 
be explained as follows: let \<f> n ) be an eigenstate of H° 
with eigenenergy e n . Then 



fx 



h.b' 



(25) 



As the Fermi energy increase the sum over e n approaches 
the sum over a complete set of eigenstates and hence fb,b' 
approaches the value (ip'bl^b')- This latter quantity van- 
ishes since the states \ifj' b ) are related to the orthogonal 
states IV'h ) by a, unitary transformation and hence re- 
main orthogonal. The lower panel of Fig. [5] shows the 
amplitude versus the switching time of the gate voltage 
for a fixed Fermi energy s-p = 0.2 a.u. and for different 
values of the applied bias. The striking feature of this 
plot is that the position of the maximum remains almost 
unchanged as function of the bias Ul- 



FIG. 7: The amplitude of the current oscillation as function 
of the switching time of the gate. The red (black) curve refers 
to the initial ground state with (without) a bound state. The 
numerical parameters are ep = 0.1 a.u., Ul = 0.15 a.u. 



As a final example, in Fig. [7] we compare the ampli- 
tude of the oscillations as function of the switching time 
t g for two different initial states with the same Fermi 
energy ep = 0.1 a.u. In one case we start, as before, 
with the constant potential Uq{x) = 0, and hence H°(x) 
does not have bound eigenstates. In the other case we 
start with a quantum well of depth Uq = —0.5 a.u. for 
|x| < 1.2 a.u.. The Hamiltonian H°\x) in this latter case 
has one bound eigenstate. A bias Ul — 0.15 a.u. in the 
left lead is suddenly switched on in both systems and 
after a time T = 150 a.u. a steady state is attained. 
For T<t<T + t g & gate voltage V g (x, t) is gradually 
switched on as before, and for t > T + t g the gate volt- 
age remains constant and equal to v g = —1.3 a.u. in the 
first case and —0.8 a.u. in the second case. Hence, both 
systems have the same asymptotic Hamiltonian H °° (x) . 
The remarkable difference between the value of the am- 
plitudes in these cases can be explained in the same way 
as in Section [III Bl 

Interestingly, in the case where the system initially has 
one bound state, the amplitude has a maximum for sud- 
den switching of the gate, i.e., t g = a.u., while in the 
case with no initial bound states the maximum appears 
at a finite value of t g . 

Similarly, we have found a maximum for small t g for 
the following situation: we start with an initial state 
without bound states. At t = a.u. we suddenly ap- 
ply a bias in the left lead and wait until a steady state is 
achieved. Then we switch on a gate in such a way that 
one bound state is created and wait until the associated 
bound-continuum transitions have decayed before we add 
another bound state to the gate with a switching time t g . 

The fact that in this case the largest amplitude for the 
current oscillations is found for switching time t g close to 
zero strongly suggests that the position of the maximum 
in the oscillation amplitude as function of t g is related 
to a transient effect. This is also supported by the fol- 
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lowing observation (see Fig[6]): the switching time t g for 
which the current oscillations are largest depends on the 
Fermi energy (for fixed bias) since the transitions from 
the bound states to the top of the Fermi sea obviously 
depend on e-p. At the same time, the position of this 
maximum is almost independent of the bias (for fixed 
Fermi energy) since the bias only leads to a slight energy 
shift for the bound states. 



IV. CONCLUSIONS 

In the theory of electron transport one usually assumes 
that the application of a dc bias to an electronic sys- 
tem attached to two macroscopic electrodes always leads 
to the evolution of a steady-state current. Recent the- 
oretical work states^ that the presence of bound states 
leads to qualitatively new features (current oscillations 
and memory effects) in the dynamics of electron trans- 
port in the long-time limit. These as well as transient 
features are investigated here in detail by numerical sim- 
ulations. In the Fourier transform of the calculated time- 
dependent current one not only finds the predicted tran- 
sitions between the bound states in the long-time limit, 
but, moreover, transitions (in the transient regime) be- 
tween the bound states and the continuum of the leads 
can also be clearly identified. We have shown that the 
amplitude of the persistent current oscillations depends 
both on the initial state and on the history of system. 
Since current and density are related via the continuity 
equation, also the time-dependent density in the long- 
time limit will therefore be history-dependent. Inter- 
estingly, these memory effects show up not only in the 
dynamical part but also in the time-independent contri- 
bution of the bound states to the density^!. 

Our results indicate that in transport calculations spe- 
cial care has to be taken if bound states are present in the 
biased system. A warning flag has already to be raised at 
the assumption of the evolution to a steady state which 
is not true in general. Of course, the theoretical anal- 
ysis predicts the existence of oscillations in the current 



but makes no statement on their relative importance as 
compared to the steady-state contribution. Our results 
show, however, that the amplitude of the oscillations lo- 
cally may very well be comparable or even larger than the 
steady-state current and therefore cannot be neglected. 
We would also like to point out that the existence of 
bound states in biased transport systems may not be an 
exotic feature in an experimental situation. For single 
molecules attached to metallic leads it is quite conceiv- 
able that some of the molecular orbitals energetically fall 
into an energy window which corresponds to an energy 
gap of the leads and those orbitals therefore cannot hy- 
bridize with any lead states and remain fully localized. 
In the case of transport experiments on quantum dots 
one could artificially create bound states by applying a 
strong attractive gate potential. 

Although our numerical simulations were performed 
for non-interacting electrons, the conclusions about the 
dynamical current oscillations apply to any effective 
single-electron theory. In particular they also apply to 
the TD Kohn-Sham equations which are in principle able 
to reproduce the time-dependent density^ (and the lon- 
gitudinal current via the continuity equation) of an in- 
teracting system if the exact exchange-correlation func- 
tional is used. Intuitively, one might expect that electron- 
electron scattering leads to a damping of the oscillations 
in the long-time limit. However, the assumption of a 
time-independent density producing a static Kohn-Sham 
potential for large times leads to a contradiction if this 
potential supports bound states since the density and 
therefore also the Kohn-Sham potential should then be- 
come time-dependent again. 
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